Explicit Scheme: Full Saint-Venant Equations

Theoretical background, numerical formulation (following K.N. Dulal's Lax Diffusive Scheme), and detailed solutions for the 1D Dynamic Wave Model.

1. Theoretical Formulation (Lax Scheme)

Following the standard procedure in K.N. Dulal's notes, the Dynamic Wave Model utilizes the full 1D Saint-Venant equations, accounting for continuity and momentum. To solve this explicitly while avoiding numerical instability, the Lax Diffusive Scheme is applied.

The governing equations (assuming no lateral inflow, $q=0$) are:

$$\textbf{Continuity:} \quad \frac{\partial A}{\partial t} + \frac{\partial Q}{\partial x} = 0$$ $$\textbf{Momentum:} \quad \frac{\partial Q}{\partial t} + \frac{\partial}{\partial x}\left(\frac{Q^2}{A}\right) + gA\frac{\partial y}{\partial x} = gA(S_0 - S_f)$$

Finite Difference Discretization

The Lax scheme replaces the spatial value at the known time level $n$ with the average of adjacent nodes to dampen numerical oscillations:

  • Time Derivative: $\frac{\partial f}{\partial t} \approx \frac{f_i^{n+1} - \frac{f_{i-1}^n + f_{i+1}^n}{2}}{\Delta t}$
  • Space Derivative: $\frac{\partial f}{\partial x} \approx \frac{f_{i+1}^n - f_{i-1}^n}{2\Delta x}$

Applying this gives the explicit formulas for Area ($A$) and Discharge ($Q$):

$$A_i^{n+1} = \frac{A_{i-1}^n + A_{i+1}^n}{2} - \frac{\Delta t}{2\Delta x} \left( Q_{i+1}^n - Q_{i-1}^n \right)$$ $$Q_i^{n+1} = \underbrace{\frac{Q_{i-1}^n + Q_{i+1}^n}{2}}_{\text{Avg}} - \underbrace{\frac{\Delta t}{2\Delta x} \left[ \left(\frac{Q^2}{A}\right)_{i+1}^n - \left(\frac{Q^2}{A}\right)_{i-1}^n \right]}_{\text{Convective Accel}} - \underbrace{g A_i^n \frac{\Delta t}{2\Delta x} \left( y_{i+1}^n - y_{i-1}^n \right)}_{\text{Pressure Gradient}} + \underbrace{g A_i^n \Delta t (S_0 - S_{fi}^n)}_{\text{Source (Bed - Friction)}}$$

The Lax Stencil

x t i-1 i (n) i+1 i, n+1

Future state at node $i$ relies on $i-1$, $i$, and $i+1$ at the current time level $n$.

Problem 1: Standard Calculation

(2078 Kartik) Rectangular channel, width = 50m, Bed slope = 0.0002, Manning’s n = 0.035, No lateral flow. $\Delta x = 1000\text{m}$ and $\Delta t = 600\text{s}$.
$Q = 38\text{m}^3\text{/s}, 36\text{m}^3\text{/s}, 35\text{m}^3\text{/s}$. $y = 1.87\text{m}, 1.82\text{m}, 1.97\text{m}$. Compute discharge and flow depth at $(i, n+1)$.

Step 1: Compute Initial Parameters at known nodes ($i-1, i, i+1$)

$$A_{i-1} = By_{i-1} = 50(1.87) = 93.5 \text{ m}^2 \quad \implies \left(\frac{Q^2}{A}\right)_{i-1} = \frac{38^2}{93.5} = 15.444$$ $$A_i = By_i = 50(1.82) = 91.0 \text{ m}^2 \quad \implies \left(\frac{Q^2}{A}\right)_i = \frac{36^2}{91.0} = 14.242$$ $$A_{i+1} = By_{i+1} = 50(1.97) = 98.5 \text{ m}^2 \quad \implies \left(\frac{Q^2}{A}\right)_{i+1} = \frac{35^2}{98.5} = 12.437$$

Step 2: Solve Continuity for Flow Depth ($y_i^{n+1}$)

$$A_i^{n+1} = \frac{93.5 + 98.5}{2} - \frac{600}{2000} (35 - 38) = 96.0 - 0.3(-3) = 96.9 \text{ m}^2$$ $$y_i^{n+1} = \frac{A_i^{n+1}}{B} = \frac{96.9}{50} = \mathbf{1.938 \text{ m}}$$

Step 3: Compute Friction Slope at node $i$

$$P_i = B + 2y_i = 50 + 2(1.82) = 53.64 \text{ m} \quad \implies R_i = \frac{91.0}{53.64} = 1.696 \text{ m}$$ $$V_i = \frac{Q_i}{A_i} = \frac{36}{91.0} = 0.3956 \text{ m/s}$$ $$S_{fi} = \frac{n^2 V_i^2}{R_i^{4/3}} = \frac{0.035^2 \times 0.3956^2}{1.696^{4/3}} = \frac{0.0001917}{2.022} = 9.48 \times 10^{-5}$$

Step 4: Solve Momentum for Discharge ($Q_i^{n+1}$)

Calculating individual terms of the momentum equation:

  • Avg Q: $\frac{38 + 35}{2} = 36.5$
  • Convective: $-\frac{600}{2000}(12.437 - 15.444) = -0.3(-3.007) = 0.902$
  • Pressure: $-9.81 \times 91.0 \times \frac{600}{2000} (1.97 - 1.87) = -892.71 \times 0.3 \times 0.10 = -26.781$
  • Source: $9.81 \times 91.0 \times 600 \times (0.0002 - 0.0000948) = 535626 \times 0.0001052 = 56.348$
$$Q_i^{n+1} = 36.5 + 0.902 - 26.781 + 56.348 = \mathbf{66.969 \text{ m}^3\text{/s}}$$

Problem 2: Wetted Perimeter Approx ($P=B$)

(2076 Bhadra) Rectangular channel, $B=25\text{m}$, $S_0=0.015$. Assume $P=B$. $\Delta x=1500\text{m}, \Delta t=600\text{s}, n=0.035$.
$Q = 30, 22, 20$. $y = 1.25, 1.45, 1.65$. Determine depth and discharge at $(i, n+1)$.
x t i-1, n Q=30, y=1.25 i, n Q=22, y=1.45 i+1, n Q=20, y=1.65 i, n+1 Q=?, y=? Δx = 1500m Δt = 10m

Step 1: Compute Initial Parameters

$$A_{i-1} = 25(1.25) = 31.25 \implies \left(Q^2/A\right)_{i-1} = 30^2 / 31.25 = 28.80$$ $$A_i = 25(1.45) = 36.25 \implies \left(Q^2/A\right)_i = 22^2 / 36.25 = 13.35$$ $$A_{i+1} = 25(1.65) = 41.25 \implies \left(Q^2/A\right)_{i+1} = 20^2 / 41.25 = 9.70$$

Step 2: Solve Continuity

$$A_i^{n+1} = \frac{31.25 + 41.25}{2} - \frac{600}{3000} (20 - 30) = 36.25 - 0.2(-10) = 38.25 \text{ m}^2$$ $$y_i^{n+1} = \frac{38.25}{25} = \mathbf{1.53 \text{ m}}$$

Step 3: Solve Momentum (Note: $P=B \implies R=y$)

$$V_i = \frac{22}{36.25} = 0.6069 \text{ m/s}, \quad R_i = y_i = 1.45 \text{ m}$$ $$S_{fi} = \frac{0.035^2 \times 0.6069^2}{1.45^{4/3}} = \frac{0.000451}{1.639} = 0.000275$$
  • Avg Q: $\frac{30 + 20}{2} = 25.0$
  • Convective: $-0.2(9.70 - 28.80) = 3.82$
  • Pressure: $-9.81 \times 36.25 \times 0.2 \times (1.65 - 1.25) = -28.45$
  • Source: $9.81 \times 36.25 \times 600 \times (0.015 - 0.000275) = 213367.5 \times 0.014725 = 3141.78$
$$Q_i^{n+1} = 25.0 + 3.82 - 28.45 + 3141.78 = \mathbf{3142.15 \text{ m}^3\text{/s}}$$ Note: The massive increase is physically consistent with the governing equations given the extremely steep slope ($1.5\%$) applied over a large time step ($10$ min) with relatively low friction.

Problem 3: Velocity-Depth Formulation

(2075 Bhadra) Wide rectangular channel ($R \approx y$). $v = 2.2, 1.8, 1.5$. $y = 1.6, 2.0, 2.4$. $S_0 = 0.01$, $n = 0.032$, $\Delta x = 1000\text{m}$, $\Delta t = 240\text{s}$. Compute velocity and flow depth at $(i, j+1)$.
x t i-1, j v=2.2, y=1.6 i, j v=1.8, y=2.0 i+1, j v=1.5, y=2.4 i, j+1 v=?, y=? Δx = 1000m Δt = 4m

Because inputs are velocities ($V$) rather than discharge ($Q$), we follow the alternate formulation (often derived directly from continuity $q = Vy$) where the momentum equation is explicitly defined in terms of $V$ and $y$:

$$\frac{\partial y}{\partial t} + \frac{\partial (Vy)}{\partial x} = 0 \quad \implies \quad y_i^{n+1} = \frac{y_{i-1} + y_{i+1}}{2} - \frac{\Delta t}{2\Delta x}(V_{i+1}y_{i+1} - V_{i-1}y_{i-1})$$ $$\frac{\partial V}{\partial t} + V\frac{\partial V}{\partial x} + g\frac{\partial y}{\partial x} = g(S_0 - S_f) \implies V_i^{n+1} = \frac{V_{i-1} + V_{i+1}}{2} - V_i \frac{\Delta t}{2\Delta x}(V_{i+1} - V_{i-1}) - g\frac{\Delta t}{2\Delta x}(y_{i+1} - y_{i-1}) + g\Delta t(S_0 - S_{fi})$$

Step 1: Solve Continuity

$$q_{i-1} = 2.2(1.6) = 3.52 \quad q_i = 1.8(2.0) = 3.60 \quad q_{i+1} = 1.5(2.4) = 3.60$$ $$y_i^{n+1} = \frac{1.6 + 2.4}{2} - \frac{240}{2000}(3.60 - 3.52) = 2.0 - 0.12(0.08) = \mathbf{1.9904 \text{ m}}$$

Step 2: Solve Momentum

$$S_{fi} = \frac{n^2 V_i^2}{y_i^{4/3}} = \frac{0.032^2 \times 1.8^2}{2.0^{4/3}} = \frac{0.003318}{2.5198} = 0.001317$$
  • Avg V: $\frac{2.2 + 1.5}{2} = 1.85$
  • Convective: $-1.8 \times \frac{240}{2000}(1.5 - 2.2) = -1.8 \times 0.12(-0.7) = 0.1512$
  • Pressure: $-9.81 \times 0.12 \times (2.4 - 1.6) = -0.9418$
  • Source: $9.81 \times 240 \times (0.01 - 0.001317) = 2354.4 \times 0.008683 = 20.443$
$$V_i^{n+1} = 1.85 + 0.1512 - 0.9418 + 20.443 = \mathbf{21.502 \text{ m/s}}$$

Problem 4: Flow Deceleration

(2073 Magh) Rectangular channel, $B=10\text{m}$, $S_0=0.0002$, $n=0.04$, $\Delta x=1000\text{m}, \Delta t=300\text{s}$.
$Q = 40, 38, 37.5$. $y = 1.9, 1.85, 2.0$. Compute discharge and depth at $(i, n+1)$.
x t i-1, n Q=40, y=1.9 i, n Q=38, y=1.85 i+1, n Q=37.5, y=2.0 i, n+1 Q=?, y=? Δx = 1000m Δt = 5m

Step 1: Compute Initial Parameters

$$A_{i-1} = 10(1.9) = 19.0 \implies \left(Q^2/A\right)_{i-1} = 40^2 / 19 = 84.21$$ $$A_i = 10(1.85) = 18.5 \implies \left(Q^2/A\right)_i = 38^2 / 18.5 = 78.05$$ $$A_{i+1} = 10(2.0) = 20.0 \implies \left(Q^2/A\right)_{i+1} = 37.5^2 / 20 = 70.31$$

Step 2: Solve Continuity

$$A_i^{n+1} = \frac{19.0 + 20.0}{2} - \frac{300}{2000} (37.5 - 40) = 19.5 - 0.15(-2.5) = 19.875 \text{ m}^2$$ $$y_i^{n+1} = \frac{19.875}{10} = \mathbf{1.9875 \text{ m}}$$

Step 3: Solve Momentum

$$R_i = \frac{18.5}{10 + 2(1.85)} = \frac{18.5}{13.7} = 1.350 \text{ m}, \quad V_i = \frac{38}{18.5} = 2.054 \text{ m/s}$$ $$S_{fi} = \frac{0.04^2 \times 2.054^2}{1.350^{4/3}} = \frac{0.00675}{1.493} = 0.00452$$
  • Avg Q: $\frac{40 + 37.5}{2} = 38.75$
  • Convective: $-0.15(70.31 - 84.21) = 2.085$
  • Pressure: $-9.81 \times 18.5 \times 0.15 \times (2.0 - 1.9) = -2.722$
  • Source: $9.81 \times 18.5 \times 300 \times (0.0002 - 0.00452) = 54445.5 \times (-0.00432) = -235.20$
$$Q_i^{n+1} = 38.75 + 2.085 - 2.722 - 235.20 = \mathbf{-197.09 \text{ m}^3\text{/s}}$$ Note: The large negative value indicates numerical flow reversal. The high roughness ($n=0.04$) and high velocity create a friction slope ($0.0045$) massively larger than the bed slope ($0.0002$), causing aggressive mathematical deceleration according to the Lax scheme.